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ABSTRACT 

We construct an analytical model of the mass distribution function of subhalos 
in a galaxy cluster in the context of the cold dark matter (CDM) theory. Our 
model takes account of two important effects; the high density of the precluster 
region and the spatial correlation of initial density fluctuations. For subhalos 
with small masses, the mass distribution function that our model predicts is 
intermediate between the Press-Schechter (PS) mass distribution function and a 
conditional mass distribution function as an extension of the PS formalism. We 
compare the results of our model with those of numerical simulations. We find 
that our model predictions are consistent with the mass and velocity distribution 
functions of subhalos in a cluster obtained by numerical simulations. We estimate 
the probability of finding large X-ray subhalos that often have "cold fronts" . The 
observed large X-ray subhalos and cold fronts may not always be the result of 
cluster mergers, but instead may be internal structures in clusters. 

Subject headings: cosmology: theory — dark matter — galaxies: clusters — galaxies: 
halos — large-scale structure of universe — X-rays: galaxies: clusters 



1. Introduction 

It is generally believed that dark matter constitutes a large fraction of the mass in the 
universe. Among various theories of dark matter, cold dark matter (CDM) theory provides 
a remarkably successful description of large-scale structure formation, and it is in good 
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agreement with a large variety of observational data. This model predicts that small objects 
are the first to form and these then amalgamate into progressively larger system. This means 
that galaxy clusters have formed via merging and accretion of galaxies or galaxy groups. 

If the dark matter is actually the main component of galaxies and clusters, iV-body 
simulations can be used to predict their bulk properties, such as position, mass, and size. 
However, it has proven difficult to find surviving substructure (subhalos) in virialized objects 
in iV-body simulations. This apparent absence of substructure, known as "the overmerging 
problem," reflects the fact that simulated subhalos are disrupted much more efficiently than 
real subhalos (Moore, Katz, & Lake 1996). 

Analytical studies showed that the overmerging problem was caused by poor force and 
mass resolutions (Moore et al. 1996). Klypin et al. (1999) estimated that the force and 
mass resolutions required for a simulated halo to survive in galaxy groups and clusters are 
extremely high: ~ 1 — 3 kpc and 10 8 — 10 9 M Q . Recently, simulations with roughly the 
required resolution have been done (Ghigna et al. 1998; Tormen, Diaferio, & Syer 1998; 
Klypin et al. 1999; Okamoto & Habe 1999; Colm et al. 1999; Moore et al. 1999; Jing & Suto 
2000; Ghigna et al. 2000; Gottlober, Klypin, & Kravtsov 2001; White, Hernquist, & Springel 
2001; Fukushige & Makino 2001; Bullock et al. 2001). These simulations certainly showed 
that subhalos survive in clusters. Okamoto & Habe (1999) compared the mass distribution 
function (MDF) of subhalos with the Press-Schechter (PS) MDF (Press & Schechter 1974) 
and a conditional MDF derived through an extension of the PS formalism (EPS; Bower 
1991; Bond et al. 1991; Lacey & Cole 1993, see §3 for detail). They found that the MDF 
of subhalos in their simulated cluster at z = is similar to the PS MDF at z — 2. It is 
also similar to the EPS MDF but the fit is a little worse. They suggested that the MDF 
of subhalos freezes out after the cluster had grown sufficiently, which in their simulation 
occurred at z ~ 2. Thus, the MDF of subhalos preserves the information about the growth 
of density fluctuations at that time. In other words, we can use observations of the MDF 
of subhalos at present to determine the hierarchical clustering which occurred in the past. 
However, the good match between the simulations and the PS MDF found by Okamoto & 
Habe (1999) is somewhat surprising as the PS MDF does not take into account the fact that 
the cluster subhalos are in a high density region. The EPS MDF does include the density 
effect, yet provides a slightly worse fit to the simulations than the PS MDF. 

Another important factor that should be taken into account when the MDF of subhalos 
in a cluster is considered is the spatial correlation of initial density fluctuations. It has been 
shown that this spatial correlation does affect the MDF of dark halos in the whole universe 
(Yano, Nagashima, & Gouda 1996; Nagashima 2001). This means that one cannot ignore 
the effect of spatial correlation on the MDF of subhalos because these subhalos originated 
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from density fluctuations that were close to one another. In this paper, we construct a model 
of the subhalo MDF in a cluster considering the effect of spatial correlations as well as the 
fact that cluster subhalos are located in a high density region. Of course, there are other 
mechanisms that could affect the MDF of subhalos especially after the cluster formation 
(e.g. tidal stripping). Although they should be taken into account in order to derive the 
MDF exactly, we do not consider them in our model for the sake of simplicity. Even so, our 
model should be useful to predict the MDF in the case where the mechanisms that occur after 
cluster formation are not very effective, and gives an upper limit to substructure if disruptive 
processes are important for a wide range of masses. Our analytic model is not affected by 
spatial and mass resolutions, which affect numerical simulations. It will be useful to compare 
with the results of numerical simulations to determine the effects of finite resolution. Even 
though some of the recent high-resolution simulations meet fully the conditions required 
to study subhalos (Klypin et al. 1999), comparing them with our model is still important. 
Such comparisons can serve as a check for the subhalo identification algorithm applied to 
the simulation results. 

Both the MDFs predicted by our model and the results from numerical simulations 
should be compared with the observations. In the past, galaxies were used to derive the 
observed MDF in a cluster. In fact, Moore et al. (1999) compared their simulation results 
with the observational data from the Virgo cluster, although they used the circular velocity 
distribution function (VDF) of galaxies instead of the MDF. They showed that the simulation 
and observed VDFs are roughly consistent with one another. However, individual galaxies 
could not be used to detect massive subhalos, because the observed circular velocity (or 
velocity dispersion) of a galaxy is at most v c ~ 300 km s -1 . On the other hand, numerical 
simulations suggest that clusters have subhalos with v c > 300 km s _1 (Moore et al. 1999). 
Until recently, it was not understood whether the lack of galaxies with v c > 300 km s _1 
requires a lack of massive subhalos in a cluster, or whether this implies that galaxies have 
lower values of v c than would characterize the massive dark matter halos surrounding the 
galaxies. Recently, X-ray observations with the Chandra Observatory have found that some 
clusters have massive subhalos moving within the clusters which have not been disrupted 
(Markevitch et al. 2000; Vikhlinin, Markevitch, & Murray 2001; Mazzotta et al. 2001, 2002; 
Sun et al. 2002). The contact surface between the cooler subhalo gas and the surrounding 
hot intracluster medium is called a "cold front" . The X-ray temperatures of these subhalos 
are ~ 4 — 7 keV, which is substantially higher that the kinetic temperatures one would derive 
from the motions of stars in the optical galaxies in the subhalos. The observations of these 
massive subhalos gives us important information about the massive end of the MDF. 

The plan of our paper is as follows. In §2, we derive the MDF in clusters analytically. 
We discuss the resulting MDFs and compare them with those from numerical simulations in 
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§3, We discuss the effect of tidal stripping and dynamical friction on the MDFs in §4. In §5, 
we derive the VDF of subhalos and compare these results with those of numerical simulations 
and observations of galaxies. We estimate the probability of finding massive X-ray subhalos 
and associated cold fronts in clusters in §6. Finally, §7 summarizes our conclusions. 



Models 



In this section, we derive the MDF of subhalos in a cluster including the effects of spatial 
correlations using the model of Yano et al. (1996). We use the so-called sharp /c-space filter 
for density fluctuations following Yano et al. (1996). We assume that the primordial density 
field is a Gaussian random field. 

First, we estimate the conditional probability, P(r, Mi, M 2 ), of finding a region of mass 
Mi with 5\ > 5 c i at a distance r from the center of an isolated, finite-sized object of mass 
M 2 , provided that the object of mass Mi is included in the object of mass M 2 (> Mi) with 
82 = ^c2 at r = 0. Here, we define 5m as the smoothed linear density fluctuation of mass scale 
M, and 5i = <5mi and 5 2 — $m 2 - Moreover, we define S c (z) as the critical density threshold 
for a spherical perturbation to collapse by the redshift z, and 5 c i = 5 c (zi) and 5 c2 = 5 c (z 2 ). 
In the Einstein-de Sitter universe, 5 c (z) = 1.69(1 + z). While Yano et al. (1996) considered 
only the case of 5 cl = 5 c2 , we consider the case of 5 c i > 5 c2 (zi > z 2 ) in general. Because the 
covariance matrix for the Gaussian distribution for two variables [<5i(r),<5 2 ] is given by 



M = 



0% 



(r) 



the probability can be written as 



P(r,Mi,M 2 ) = p[5i>5 c i\5 2 = 5 c2 ] 

= p[Si > 6 cl , 5 2 = 5 c2 ]/p(5 2 = 5 c2 ) 

1 f 00 f-h-e(r> 2c ] 2 
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where v\ and v 2 are defined by 
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respectively (Yano et al. 1996). In equations (3), a(M) is the rms density fluctuation 
smoothed over a region of mass M, and is given by 
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In this equation, V is the volume, k is the wave number, and 5k is the Fourier components 
of density fluctuations. The critical wave number, k c (M) is the wave number corresponding 
to mass scale M and is given by k c (M) = 7r/(3M/47rp), where p is the current mean mass 
density of the universe. Moreover, e(r) is defined by 



e(r) = 



^(r) = <(J2(ro)«5i(r + r)> = 



y fkc{M 2 ) 
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Jo 



2 sm(kr ) A _ ]2 

! kr 
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(5) 
(6) 



which corresponds to the two-point correlation function. 
We can rewrite equation (2) as 

1 



where 
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y/l -e 2 (r) y/1 - e 2 (0)a 2 (r) <?i 
and a(r) = e(r)/e(0). The spatially averaged conditional probability for r < i? is defined as 
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The definition of P(M 1 , M 2 ) (equations [7] and [9]) is similar to the probability that the 
density contrast on a scale of M exceeds S c in field: 



1 f°° 

f(>6 c ,M) = — 

V Zn Js c /a 

From this, the PS MDF can be obtained as 
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h PS dM = 2 



M 



df 



dM 



dM , 



(11) 



(Press & Schechter 1974). 



Using the similarity between /(> <5 C , M) and P(r, M 1; M 2 ), we define the MDF of the 
collapsed objects (subhalos) in the region of a mass scale of M by differentiating P(M, M ) 
and multiplying it by 2 (PS approximation): 



n S p S dM = 2 



Mo 
M 



dP(M, M ) 
dM 



dM 



(12) 
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This expression is based on the usual PS assumptions. However, we emphasize that the 
MDF (equation [12]) includes the effect of the spatial correlation of fluctuations through the 
two-point correlation function (equation [6]). Moreover, it is the MDF in a high density 
region, rather than the MDF for an average region of the universe. As will shown in the 
next section, our model is a further extension of the EPS MDF by taking account of spatial 
correlations of halos. Thus, we refer to this MDF as SPS MDF. We assume that the region 
of the mass scale of M will collapse into a cluster and the subhalos in the region also fall 
into the cluster. The density fluctuations in the cluster region should stop growing as the 
cluster grows, because the kinetic energy released at the cluster collapse allows the subhalos 
to move in the cluster, which prevents their further growth. We call the redshift when this 
occurs the effective formation redshift of the cluster, Zf. Thus, we use 5 C \ — S c (zf) or z\ — Zf 
in equations (2) and (12). 

Strictly speaking, the conditional probability P(M 1 ,M 2 ) derived in this section is not 
sufficient to obtain the exact MDF. This is because an isolated object of mass scale M 2 
must have the maximum peak density with 5 2 = S c2 at r = 0. Thus, we need an additional 
condition. However, Yano et al. (1996) showed that the effect is negligible as long as we 
consider a fluctuation power spectrum with an index of n < —2 and the mass scale with 
a(M) > 1. Our calculations in the following sections satisfy these conditions. 



3. Mass Distribution Function 

We calculate the MDF of subhalos in a cluster, nsps, for a CDM universe with cos- 
mological parameters of Qq = 1, A = 0, <7g = 0.7, T = 0.5, and h = 0.5, where a 8 is the 
rms mass fluctuation on 8/i" 1 Mpc, and T is the so-called shape parameter of the CDM 
spectrum. We use a Hubble constant of Hq = lOO/i km s" 1 Mpc" 1 . The cluster mass 
at z = is M = 4.3 x 10 14 M Q . These parameters are the same as those used for an 
ultra-high resolution numerical simulation done by Ghigna et al. (2000). We calculate the 
conditional probability, P(M, M ), spatially averaged in a precluster region, r < R = R , 
where R = (SMo/Airp) 1 ^ 3 . From now on, we assume R = R unless otherwise stated. 

For comparison, we also calculated a conditional MDF, that is, the number of halos 
with mass between M and M + dM at Zf that are in a halo of mass M at z (M < M , 
z < zj). It is given by 



n EPS (M, z f \M , z )dM = M±f(S,u\S , u ) 



dS 



dM 



dM , (13) 
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where 



f(S,tu\S ,u )dS = 



UJ — Co>o 



exp 



(u - cup) 2 
'2(S-S )_ 



dS , 



(14) 



(27T)V2( > S_ > S )3/2 

5 = a 2 (M), 5 = (T 2 (M ), u = Sc(zf), and cu = 8(z ) (Bower 1991; Bond et al. 1991; Lacey 

6 Cole 1993). This formalism is often called the extended Press-Schechter formalism (EPS). 

We note that the difference between nsps and nsps comes from whether the effect of 
the spatial correlation of fluctuations is included or not. If we ignore the spatial correlation, 
P(Mi, M 2 ) should be written as 



P(M 1 ,M 2 ) = p(8 1 >8 cl \8 2 = 8 c2 ) 
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(see Yano et al. 1996). In this case, it can easily be shown that the right hand of equation 
(12) is identical to that of equation (13) (see Bower 1991). We also note that at r = 0, the 
right hand of equation (2) is identical to that of equation (15) (Yano et al. 1996). We refer 
the right hand of equation (15) as P E ps(^i, M 2 ), so that 



n EPS dM = 2^ 



dP EPS (M,M ) 



dM 



dM . 



(17) 



We would like to point out that the equation (9) gives the probability distribution of the 
mass of halo progenitors as Peps(^1; M 2 ). Thus, it could be possible to construct dark halo 
merger trees including the spatial correlation of dark halos, which would be important for 
semi-analytic models of galaxy formation. 

We also compare nsps with the conventional PS MDF (equation [11]). Instead of n PS , 
we use the MDF, 



n PS dM = —hp S dM = 2— 
p M 



df 



dM 



dM . 



(18) 



Of course, the PS MDF includes neither the effect of the spatial correlation of fluctuations 
nor that of the density excess in a precluster region. 



Figure 1 shows the three MDFs at z = z 

1 dN 
VdM ~~ 



■ 0. For the vertical axes, we use 
V ' 



(19) 



- 8- 



where V = AttR 3 c /3 and i = SPS, EPS, and PS following Ghigna et al. (2000). We use 
capital N for cumulative numbers from now on. The cluster radius, R c , is the same for all 
the three models and is given by R c = (3M /47rp ) 1//3 , where p is the average density of the 
cluster and is 178 times the critical density of the universe for our cosmological parameters. 
In equation (10), we use 4 = S c (zf). 

As the effective formation redshifts of a cluster, we use Zf = 0, 0.5, and 1 for nsps 
and nps- On the other hand, if we chose Zf = in the EPS MDF, we would find trivially 
that ngps = except at M = M . (That is, there are no subhalos except for the cluster 
itself.) Thus, for the EPS MDF we show the results for Zf = 0.1, 0.5, and 1. The mass 
distribution function nsps increases with Zf for smaller M, which is also true of tzeps and 
n PS . Figure 1 also shows that n EPS < n SPS < n PS for Zf ~ = z and M < 10 12 M G . This 
can be understood as follows: When Zj- C 1 (or 5 c [zj] 5 c [z ]), M <^ M (or e[0] <^ 1), and 
6 c (z f )/cr(M) < 1, we obtain 



and 



This implies that 



a(M) 
J a(M) 



1 - 



4(0) 
4(0) 



a(r) 



$c(Zf) 



(3'<(3< 5 c (z f )/a(M) < 1 , 
from equations (8) and (16) because < a(r) < 1. Thus, we obtain 
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dM a(M) 
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from equations (7), (15), and (10). Since the relation 



dP EPS (M,M ) 



dM 
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(22) 

(23) 

(24) 
(25) 
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always holds for < r < R, we have 



^eps < n SPS < n PS 



(27) 
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from equations (12). (17), and (18). Relation (27) reflects the fact that in the model of Yano 
et al. (1996) the critical density threshold of the collapse of the cluster is effectively between 
zero and 5 C (0) at < r < R. On the other hand, when heps is derived, it is implicitly 
assumed that the critical density threshold is uniform and is S c (0) in the precluster region. 
These results show that at least ueps and nps are not appropriate to describe the MDF of 
subhalos especially when Zf < 1. 

On the other hand, when z f > 1 (or S c [z f ] > 5 c [z ] = 1.69), M < M (or e[0] < 1), we 
obtain 

P(r) « (3' « 5 c (z f )/a(M) , (28) 

from equations (8) and (16). This means that npps ~ ^sps ~ ^ps- Figure 1 shows that this 
relation is satisfied even when Zf ~ 1 for M < 10 12 M Q . 

We compare the above results with those of ultra-high resolution numerical simulations 
done by Ghigna et al. (2000). In their simulation, a cluster contains ~ 5 million particles 
within the final virial radius and is simulated using a force resolution of 1.0 kpc (0.05% of the 
virial radius). The particle mass is 1.1 x 10 8 M & . Of the three MDFs nsps, ^eps, and nps, 
only n S ps includes both the effects of the spatial correlation of fluctuations and the density 
excess in a precluster region, so we concentrate on nsps- Figure la shows that the result of 
numerical simulation at z = is consistent with nsps at Zf ~ 0.5. The change of the MDF 
slope of the simulated subhalos at M < 10 9 M seems to be caused by incompleteness due to 
limited resolution. The formation redshift, zj, is reasonable because the cluster in Ghigna et 
al. (2000) does not undergo major mergers since z = 0.5 and the virial mass increases only by 
< 30% since that epoch. In fact, Ghigna et al. (2000) compared the MDF at z = with that 
of subhalos at z — 0.5 within the same physical volume, as determined by the virial radius 
of the cluster at z — 0, and found that they are almost identical. Thus, the MDF is "frozen 
in time" since z = 0.5, and Zf should be > 0.5. Note that since the comparison was done 
within the same physical volume, the accretion of halos from outside cluster does not affect 
the evolution of the MDF. Simulation of clusters that collapse more recently (zf ~ 0) would 
be useful to confirm whether nsps is actually superior to npps and nps, because the three 
MDFs are more disparate at smaller Zf (Figure 1). Moreover, a recently formed cluster, the 
effects we have not considered (e.g., tidal stripping in the cluster) would not have much time 
to operate. One can estimate the formation redshift Zf of the simulated cluster by finding 
the redshift at which the MDF (or VDF; see next section) stop evolving. 

We also compare our model with another numerical simulation done by Okamoto & 
Habe (1999). They adopted cosmological parameters of Q = 1, A = 0, er 8 = 2/3, T = 0.5, 
and h = 0.5. The present mass of a cluster is M = 9.3 x 10 14 M . In their simulation, the 
cluster contains one million particles within the final virial radius and is simulated using a 
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force resolution of 5.0 kpc (0.2% of the virial radius). The particle mass is 1.08 x 10 9 M Q . 
Thus, the resolution is a little worse than Ghigna et al. (2000). Figure 2 shows the MDF 
of subhalos in the cluster. For vertical axes, we use dN/dM = (i =SPS, EPS, and PS) 
following Okamoto & Habe (1999). For Zf > 1, the three models are not very different for 
M < 10 12 M , which is the same as in Figure 1. The result of Okamoto & Habe (1999) at 
z = is also shown in Figure 2. As can be seen, nsps at Zf ~ 2 is consistent with the MDF 
of subhalos obtained by Okamoto & Habe (1999). In fact, they indicated that the MDF 
does not evolve significantly from z = 2 to the present; this cluster seems to have formed 
earlier than that in Ghigna et al. (2000). Okamoto & Habe argue that the reason why the 
MDF in their simulated cluster agrees better with n PS than with ueps at Zf ~ 2 may be 
because of their halo-finding algorithm. However, our results suggest this difference may be 
because tieps does not include the effect of spatial correlations of fluctuations and because 
nps is almost the same as nsps at Zf ~ 2 (Figure 2). 



4. Tidal Stripping and Dynamical Friction 

We should note that we do not consider some mechanisms that could change the MDF of 
subhalos in a cluster after the cluster formation; these mechanisms might have been effective 
even in cluster progenitors. 

First, tidal stripping reduces the mass of subhalos. We assume that the density profiles 
of subhalos and the main cluster are represented by the so-called NFW profile: 

P(r) = ? P f , 2 , (29) 
r(r + r s y 

where r s and p are the characteristic radius and density of the halo, respectively (Navarro, 
Frenk, & White 1997). The mass profile is given by 

M(r) = M vi J(x)/f(C) , (30) 

where 



and 



/(x)=ln(l + x)--^-, * = -, (31) 
1 + x r s 

r ■ ( M ■ \ - a084 

C=^124 '7" (32) 
r s \lh- l M & ) 

(Klypin et al. 1999). In this section, we often refer to subhalos and the main cluster as 
halos, and represent the total mass of halos by their virial masses M vir ; their virial radii are 
represented by r vir . The shape of the profile is characterized by the concentration parameter 
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C. The circular velocity of a halo at a radius r is defined as v c (r) = [GM(< r)/?"] 1 / 2 . For 
the NFW profile, the maximum of the circular velocity occurs at r max 2r s . and it is given 
by 

"GM vir /(2) 1 1/2 



v 



peak 



2/(C) 



(33) 



The tidal radius r t of a subhalo with mass M(r) and peak circular velocity t> pca k moving 
at a radius i? from the center of a cluster with mass M (R) and peak circular velocity v is 
given by the radius at which the gravity force of the subhalo is equal to the tidal force of 
the main cluster. The radius can be derived by solving the following equation for the tidal 
radius r t : 

R\ 3 M(r t ) o R dM 



r t J M (R) M OR 



which is equivalent to 

f(x t ) ( x t V / r s v 



f(x R ) \x R J \R s v peak 



2 - 



x 



R 



l + x R ff{x R ) 



(35) 



where x t = r t /r s , xr = R/R s , and R s is scale radius of the cluster (Klypin et al. 1999). 

On the other hand, tidal stripping can occur by resonances between the force the subhalo 
exerts on the dark matter particle and the tidal force by the main cluster. In this case, the 
tidal radius can be obtained by solving the following equation for the tidal radius: 

GM(r t ) _ GMq(R) 



which is equivalent to 



f(x t ) = ( xt\ 3 / r s v \ 
f(x R ) \Xr) \R s v pe3k J 



(37) 

(Klypin et al. 1999). We take the smaller of the two estimates of r t . 



Figure 3 shows the mass fraction of a subhalo within the tidal radius as a function of the 
'subhalo' virial mass M vir . The cosmological and cluster parameters are the same as those 
in Figure 1. As can be seen, tidal stripping is very effective when the distance from the 
cluster center is small. Moreover, it is effective when M vir is large, because the concentration 
parameter C is small (equation [32]). In a real cluster, the effect of tidal stripping depends 
on the orbit of a subhalo in a cluster. Since it is difficult to discuss the orbits of subhalos 
analytically, numerical simulations are required to estimate the effect of tidal stripping on 
entire subhalos in a cluster quantitatively. Ghigna et al. (1998) estimated that for a circular 
velocity larger than ~ 10 2,2 km s -1 (M vir ~ 2.5 x 10 11 M ), the masses of subhalos are 
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smaller than those of isolated halos by ~ 20 — 50% on the average due to tidal stripping. 
Thus, the tidal stripping would affect the comparison between our numerical model and the 
results of numerical simulations done in §3, especially for the high mass end of the MDFs. 

Second, the most massive subhalos often merge with the central "cD galaxy" via dy- 
namical friction. Thus, we may overestimate the number of these subhalos. The time-scale 
of dynamical friction can be estimated by using Chandrasekhar's formula. When the density 
distribution of a cluster is approximated by an isothermal density distribution, a subhalo 
with the initial position of R reaches the cluster center in the time of 

t 117 R2v ° 



M(r t ) 1 



10 13 M, 



© 



(39) 



where A = R c M (R)/[RM(t t )], and R c is the cluster radius (Binney & Tremaine 1987). Note 
that the time-scale (equation [39]) is not much different from the one obtained by Klypin et 
al. (1999) without the approximation of an isothermal distribution. 

Figure 4 shows the time-scale of dynamical friction £ fric as a function of the subhalo virial 
mass. The cosmological and cluster parameters are the same as those in Figure 1. Massive 
subhalos near the cluster center are most affected by dynamical friction. The fraction of 
subhalos that fall into the cluster center depends on the distribution of their orbits. In 
the numerical simulations done by Ghigna et al. (2000), 60% of subhalos with the circular 
velocity of > 300 km s" 1 (M vir > 2.5 x 10 12 M ) have merged into the central cD galaxy 
since z — 3, although the mergers have been less effective after z — 1. Thus, the MDFs 
derived in §3 should be compared with the results of numerical simulations and observations 
with caution for those massive subhalos. 

On the other hand, the distribution of peak circular velocities v pea k is less affected by 
tidal stripping, because the removal of the outer parts of a subhalo does not affect very 
strongly. Klypin et al. (1999) showed that even if ~ 70% of the initial mass of a subhalo is 
removed by tidal stripping, t> pca k changes only by ~ 20%. Since the numerical simulations 
done by Ghigna et al. (1998) showed that tidal stripping reduces the mass of massive subhalos 
at most by ~ 50%, the tidal stripping may not affect the distribution function of values of 
Vpcak (VDF) for subhalos significantly. Thus, in the next section, we investigate the VDF 
for subhalos, although the high-velocity end of the VDF may still be affected by dynamical 
friction, or mergers of subhalos. Observationally, w pea k can be obtained more easily than M 
through the velocity dispersion of stars or the rotation curve of a galaxy. Thus, the VDF is 
more useful for comparison with observations than the MDF. 
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5. Velocity Distribution Function 



In order to derive the VDF of subhalos from the MDF, we need the relation between the 
velocity and mass. Since there is no simple analytically derived relation, we use the relation 
derived by numerical calculations for isolated halos. Navarro, Frenk, & White (1997) found 
that the relation between w pea k and M is 



for the Einstein de-Sitter universe. Although the relation is for z — 0, it is not much different 
for z < 2 (Ghigna et al. 1998; Bullock et al. 2001). We assume that this relation holds for 
subhalos, at least before the collapse of the cluster at z = Zf. Using equation (40), we 
transform the MDFs derived in §3 into VDFs; the results are shown in Figure 5. Note that 
when we discuss the VDFs, n» (i = SPS, EPS, and PS) is a function of t> pca k, and 



Since f pea k is a monotonically increasing function of M, the relation among nsps, ^eps> and 
nps is the same as that for MDF. In Figure 5, the results of numerical simulations are also 
shown (Ghigna et al. 2000). In general, nsps at Zf = 0.5 is consistent the result of numerical 
calculation as was shown to be the case for the MDF. This may mean that tidal stripping 
does not significantly affect the MDF of subhalos in the simulated cluster. The deviation 
of the simulated cluster data from the model at w pea k < 10 L9 km s _1 is caused by finite 
numerical resolution. The slope of the VDF in the simulated cluster is a little steeper than 
that of n SPS for f pea k > 10 2 ' 3 km s" 1 if we remove the largest w pea k point. This may be caused 
by mergers of massive subhalos with the central cD galaxy. 

If the VDF is not significantly influenced by tidal stripping, we can use it to derive the 
effective formation redshift of a cluster. Figure 6 shows the observed cumulative number 
of subhalos in the Virgo cluster derived from galaxy data (Binggeli et al. 1985; Moore et 
al. 1999) and the integrated nsps for various Zf. The cosmological parameters for nsps are 
the same as those in Ghigna et al. (2000). The comparison suggests that Zf ~ 2, although 
the slope of the Virgo data does not match that of nsps for ^peak/^o ^0.1. Here, vo is the 
peak velocity t> pca k of the rotation curve for applicable to the total mass of the cluster. The 
deviation for low f pea k subhalos may be because the corresponding galaxies are too faint to 
be detected. Note that that the derived formation redshift, Zf ~ 2, is the redshift at which 
the growth of subhalos stopped. This corresponds to the time when the potential wells of 
the cluster progenitors, in which the subhalos were located, became sufficiently deep, even 
if the progenitors had not yet merged into one cluster. Since the subhalos stopped growing, 




(40) 
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the formation of galaxies in them is also expected to decrease for z < Zf. In the irregular 
Virgo cluster, Zf may correspond to the time at which the oldest galaxies were formed, even 
though the cluster itself is still merging at present. 

We note that we cannot rule out the possibility that Zf is somewhat smaller than 2 
because of the ambiguity of the number and circular velocity of massive subhalos. In general, 
for massive subhalos, dynamical friction may cause them to merge into one large (cD) galaxy, 
which reduces their number. However, although there is a moderately large galaxy (M87) at 
the center of the Virgo cluster, it is not usually classified as a cD galaxy and is probably not 
the most luminous galaxy in the cluster. Alternatively, the circular velocity of the galaxy 
residing in a massive subhalo may not represent that of the massive subhalo; the former 
may be smaller than the latter. In fact, Davis & White (1996) showed that the temperature 
of diffuse X-ray gas is systematically higher than that of kinetic temperature of stars for 
bright elliptical galaxies. If the temperature of the X-ray gas is considered to represent the 
depth of potential well of the halo around the elliptical galaxy, it shows the inconsistency of 
circular velocity between galaxies and subhalos. Moreover, Matsushita (2001) showed that 
bright elliptical galaxies with extend X-ray emission tend to have high X-ray temperature 
compared to those with compact X-ray emission for a given stellar velocity dispersion. On 
the other hand, the X-ray gas may be affected by other heating processes (e.g., supernovae). 

These results suggest the existence of very large, massive subhalos whose X-ray temper- 
ature is much higher than kinetic temperature of the stars in the galaxies located in them. 
Recent Chandra observations have revealed the presence of such massive subhalos in clus- 
ters; most of them are moving with the Mach number of ~ 1 in clusters without losing their 
identity. Thus, in the next section, we discuss the relation between the large X-ray subhalos 
and VDF. 



6. Cold Fronts 

Large X-Ray subhalos are found in some clusters as low temperature X-ray components 
with cold fronts (Markevitch et al. 2000; Vikhlinin et al. 2001; Mazzotta et al. 2001, 2002; 
Sun et al. 2002). The cold fronts are contact interfaces between the cooler gas in the subhalo 
and the hotter intracluster medium in the main cluster. The ratio of X-ray temperature of 
the subhalos to that of their host clusters is about 1/2. Thus, if the X-ray temperature is 
proportional to fp Cak , the corresponding velocity ratio is Vpe^/vo ~ l/y/2. 

Mazzotta et al. (2001) suggested that a cluster with a large subhalo may be the result 
of the collapse of two different perturbations in the primordial density field on two different 
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linear scales at nearly the same location in space. As the density field evolves, both per- 
turbations start to collapse. The small-scale perturbation collapses first and forms a large 
subhalo, and then the larger perturbation collapses and forms a cluster. If the initial posi- 
tion of the subhalo was only slightly offset from the center of the cluster, it would infall into 
the cluster center with small velocity and oscillate around the cluster center. This would 
prevent the cool gas in the subhalo from being removed by ram-pressure stripping or the 
Kelvin-Helmholtz instability. We note that Fujita et al. (2002) showed that the observed 
subhalos with cold fronts are stable for 0.5 — 1 Gyr against large-scale Kelvin-Helmholtz 
instability. 

Therefore, it would be interesting to estimate the probability of finding such large sub- 
halos near centers of clusters, and to compare this probability with the observed rate of 
occurrence of cold fronts. If we had enough observational data about the massive subhalos, 
these could be used to study the validity of our MDF or VDF model. However, at present the 
number of clusters in which large X-ray subhalos are found is too small to allow a detailed 
discussion of the rate of occurrence of cold fronts. Thus, in this study, we just predict the 
nature of the massive subhalos, assuming that our model is correct. The cumulative number 
of the most massive halos is larger for smaller Zf (Figure 6), because hierarchical clustering 
proceeds furthest. Thus, we use nsps at Zf = to estimate the maximum probability of 
finding those subhalos. 

Figure 6 shows that the cumulative number, iV(0.71 < fpcak/^o < 1)> is about 0.3, which 
means that at most one third of clusters should have large subhalos, some of which may have 
cold fronts similar to the observed ones. In order to find the initial position of the subhalos 
in a precluster region, we calculate the number by changing the radius, R, within which the 
conditional probability (equation [9]) is averaged. Figure 7 shows the result for the same 
cluster in Figure 6. Since N(r < R, 0.71 < v pcak /v < 1) is almost constant for R/R > 0.5, 
it is most likely that the initial position of the subhalo is -R/-R0 ^ 0.5. Because we do not 
consider the spatial distribution of mass in a precluster region and the detailed evolution of 
cluster collapse, we cannot exactly predict the infall or oscillating velocity of the subhalo in a 
cluster. However, if the initial position of the subhalo is well inside the precluster region, the 
infall or oscillating velocity may be subsonic or at most transonic, which may be consistent 
with the observations (Markevitch et al. 2000; Vikhlinin et al. 2001; Mazzotta et al. 2001, 
2002). Note that if the initial position of a subhalo is near r = R , the Mach number of the 
infall velocity should be ~ 2 — 3 (Sarazin 2002). The above result suggests that most of the 
X-ray substructures in clusters have not fallen in from large distances, as would be true for 
cluster mergers. Instead, most of the cold fronts may be due to internal subhalos. The fact 
that clusters having subhalos moving with highly supersonic velocities (e.g. 1E0657— 56; 
Markevitch et al. 2002) are rarely found may support this idea, although the number of 



-16- 



clusters for which the velocities of subhalos are derived is small. 

As is mentioned in §4, the large subhalos should be affected by dynamical friction 
(v c ^ 300 km s" 1 according to Ghigna et al. 1998). In approximately the time shown in 
equation (39) the subhalos reaches the cluster center. If the subhalos are not completely 
disrupted before they reach the cluster center, the clusters should have a potential well with 
two different spatial scales. In fact, Ikebe et al. (1996) found a double dark halo distribution 
in the Fornax cluster. On the other hand, in the central regions of many clusters, cool gas 
is observed (Fabian 1994). These subhalos may provide cool gas into the cluster center. 

7. Conclusions 

We have constructed a model of the mass distribution function (MDF) of subhalos in 
clusters, including the fact that the subhalos grew from spatially correlated density fluctua- 
tions in high-density regions. If a cluster forms at a low redshift, the derived mass distribution 
function of subhalos is intermediate between the Press-Schechter mass distribution function 
and a conditional mass distribution function derived by an extension of the PS formalism. 
However, if a cluster forms at high redshift, our subhalo MDF is nearly identical to these 
other MDFs for low mass subhalos. The velocity distribution functions show similar trends. 
We compare our model MDF with the results of iV-body numerical simulations. We find 
that they are consistent with each other. The comparison of our model with the subhalo 
population derived from galaxy data suggests that the subhalos in Virgo cluster stopped 
growing at z ~ 2. Our model predicts that at most 1/3 of clusters should have large X-ray 
subhalos like the ones observed as cold fronts with Chandra. Many of the observed X-ray 
substructures in clusters may not be the result of cluster mergers; instead, they may be due 
to internal substructures within the clusters. 
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and T. Okamoto for useful comments. We acknowledge the referee for helpful suggestions. 
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Fig. 1. — Subhalo mass distribution functions for Zf = or 0.1 (dotted lines), Zf = 0.5 (solid 
lines), and Zf = 1.0 (dashed lines), (a) nsps, (b) ueps, and (c) nps- The MDF at z = 
obtained from the numerical simulation by Ghigna et al. (2000) is shown by the circles. 
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Fig. 2. — Subhalo mass distribution functions for Zf = or 0.1 (dotted lines), Zf = 1.0 (solid 
lines), and Zf = 2.0 (dashed lines), (a) nsps, (b) heps, and (c) nps- The MDF at z = 
obtained from the numerical simulation by Okamoto & Habe (1999) is shown by the circles. 
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Fig. 3. — The fraction of mass within the tidal radius of subhalos with different masses for 
several distances from the cluster center (R). The cosmological and cluster parameters are 
the same as those in Figure 1. 




log Kir /H> 

Fig. 4. — The dynamical friction time for different masses of subhalos for several distances 
from the cluster center (R). The cosmological and cluster parameters are the same as those 
in Figure 1. 
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Fig. 5. — Velocity distribution functions for Zf = or 0.1 (dotted lines), Zf = 0.5 (solid 
lines), and Zf = 1.0 (dashed lines), (a) n SPS , (b) n E psj and (c) n PS . The VDF at z = 
obtained from the numerical simulation by Ghigna et al. (2000) is shown by the circles. 
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Fig. 6. — Circles show the observed cumulative number of subhalos in the Virgo cluster 
derived from galaxy observations (Binggeli et al. 1985; Moore et al. 1999). The lines give 
the integrated SPS MDF for z f = or 0.1 (dotted lines), z f = 1.0 (solid lines), and z f = 2.0 
(dashed lines). The integrated SPS MDF was calculated using the cosmological parameters 
in Ghigna et al. (2000). 
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Fig. 7. — The average number of subhalos with 0.71 < v pcak /v < 1 within r < R. We 
assume that Zf — 0. 



